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Abstract 

Procedures  are  presented  for  the  systematic  design  of  digital 
filters  that  contain  poles  and  zeros.  The  procedures  are  simple,  fast, 
and  effective.  All  of  the  important  algorithms  are  of  the  Levinson-type. 
The  first  key  idea  in  the  paper  is  that  one  may  begin  a design  by  posing 
a linear  prediction  problem  for  a stochastic  sequence.  The  second  is 
that  a high-order  '^whitening^  filter  may  be  constructed  for  this  sequence 
and  '^Inverted^  to  yield  a high-order  all-pole  filter  whose  spectrum 
approximates  the  spectrum  of  the  stochastic  sequence.  The  third  key 
idea  is  that  the  all-pole  filter  may  be  used  to  generate  consistent 
unit  pulse  and  covariance  sequences  for  use  in  the  Mullis-Roberts 
algorithm.  This  algorithm  is  then  used  to  obtain  a low-order  digital 
filter,  with  poles  and  zeros,  that  approximates  the  high-order  all-pole 
filter.  The  results  demonstrate  that  the  Mullis-Roberts  algorithm, 
together  with  the  design  philosophy  of  this  paper,  may  be  used  with 
profit  to  reduce  filter  complexity  and  to  design  spectrum-matching 
digital  filters. 
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1.  Introduction 

This  is  a paper  about  the  systematic  design  of  autoregressive-moving 
average  (ARMA)  digital  filters.  There  are  three  key  ideas  in  the  paper. 

The  first  is  that  one  may  begin  the  design  of  an  ARMA  digital  filter  by 

specifying  the  power  spectral  density  of  a stochastic  sequence  and  then  j 

1 

posing  a classical  linear  prediction  problem  for  the  stochastic  sequence.  | 

j 

The  second  is  that  this  prediction  problem  may  be  solved  by  designing  a \ 

high-order  moving  average  (MA)  prediction  filter.  This  prediction  filter  ' 

is  related  to  an  MA  whitening  filter  which  may,  in  turn,  be  "inverted"  to 

give  an  autoregressive  (AR)  filter  termed  an  inverse  filter.  Identical 

argumentation  lies  at  the  heart  of  AR  or, maximum  entropy  (ME)  spectrum 

analysis  and  explains  why  the  magnitude-squared  frequency  response  of  the  ^ 

AR  filter  should  approximate  the  power  spectrum  of  the  stochastic  sequence 

[1].  The  appropriate  design  equations  are  linear  equations  (termed  normal 

equations)  that  may  be  efficiently  solved  using  Levinson-type  algorithms. 

For  our  purposes  the  value  of  the  high-order  AR  filter  is  simply  that  it 

provides  a handy  mechanism  for  generating  consistent  unit  pulse  and  covari- 
00  00 

ance  sequences  ®nd  respectively.  These  sequences  approximate 

the  unit  pulse  and  covariance  sequences  of  an  idealized,  unrealizable, 
digital  filter. 

00  00  j 

The  third  key  idea  is  that  the  sequences  and  „ ™ay  be  used  l 

in  the  approximation  algorithm  published  by  Mullls  and  Roberts  in  their 

remarkable  paper  on  the  use  of  first-  and  second-order  information  in  i 

00  00  ^ 
discrete-time  system  design  [2],  From  and  for  the  high-order  j 

M N 

AR,  the  finite-length  sequences  and  are  used  in  the  Mullis-  1 

Roberts  algorithm  to  design  an  ARMA  (N,M)  digital  filter.  The  notation  { 
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ingly,  AR  (N^^)  will  denote  an  all-pole  digital  filter  with  poles  and 
MA  (M^)  will  denote  an  all-zero  digital  filter  with  zeros. 

Our  results  indicate  that  the  design  procedure  outlined  above  is 
simple,  effective,  and  fast.  All  of  the  important  algorithms  are  Levinson- 
type  algorithms.  Moderate-order  filters  (e.g.,  N=16  and  M=16)  may  be  de- 
signed in  approximately  8 seconds  of  CDC  6400  CPU  time.  Design  parameters 
are  center  frequency,  bandwidth,  and  stopband  rejection.  Even  for  low- 
order  filters  (e.g.  N=8  and  M=8)  one  may  achieve  60  dB  rejection  over  a 
transition  band  that  is  a small  percentage  of  the  foldover  frequency  with 
approximately  3dB  passband  ripple. 

As  an  approximation  tool  the  Mullis-Roberts  algorithm  is  very  effec- 
tive: for  example,  an  ARMA  (N,M)  filter  of  complexity  C = N + M is  shown  in 

many  instances  to  be  a very  good  approximation  to  an  AR  (N^^)  filter  for 
C <<  Designs  for  = 256  and  C = 32,  and  N = 64  and  C = 16  illus- 

trate the  point.  Our  experience  Indicates  that  the  designs  presented 

here  are  superior  to  those  achieved  by  simply  matching  the  impulse  sequence 
1W“M 

{h^}^  using  the  Burrus-Parks  algorithm  [3].  Furthermore  there  is  no 
difficulty  with  stability  because  both  the  AR  (N^^)  and  ARMA  (N,M)  filters 
are  guaranteed  to  be  stable. 

Preliminary  design  efforts  along  the  lines  of  this  paper  were  reported 
in  [4].  The  designs  were  begun  with  a moderate-order  MA  (M^^)  filter  (e.g., 
Mj^  = 32)  designed  using  the  methodology  of  [5].  The  MA  filter  was  used 
to  generate  consistent  unit-pulse  and  covariance  sequences  for  use  in  the 
Mullis-Roberts  algorithm.  The  results  reported  were  mixed  because  the 
unit-pulse  and  covariance  sequences  were  not  well-suited  for  subsequent 
approximation  using  the  Mullis-Roberts  algorithm.  In  our  opinion,  and 
the  results  seem  to  corroborate  this  view,  one  should  begin  ARMA  (N,M) 


approximations  with  very  high-order  AR  (N^^)  designs  (e.g.,  = 256) 

rather  than  very  long  or  short  MA  designs.  The  reasons  are  three-fold; 

(!)  AR  power  spectrum  approximatrons  are  well-understood  [6]  and 

virtually  free  in  terms  of  numerical  design  effort,  (ii)  the  sequences 
00 

and  {r,  } are  obtained  from  simple  regression  equations,  and  (iii) 

R -oo  ^ 

the  AR  sequences  and  are  infinite-length  and,  therefore,  appar- 

ently better  suited  to  subsequent  matching  and  approximation  by  Infinite 
length  ARMA  unit-pulse  and  covariance  sequences. 

We  remark  that  the  approximation  algorithm  of  Mullis  and  Roberts 
gives  exact  matching  of  the  AEMA  (N,M)  and  AR  (N^^)  unit  pulse  sequences 
over  M + 1 indices  and  only  approximate  matching  of  the  covariance 
sequences  over  N + 1 indices.  When  oiS.y  a covariance  sequence  is  to  be 
matched,  then  the  Mullis-Roberts  algorithm  reduces  to  the  Levinson  algo- 
rithm for  obtaining  an  exact  covariance  match  over  N + 1 indices. 


We  feel  the  statistically-related  filter  design  procedures  reported 
here  and  in  [5]  and  [7]  are  becoming  well-enough  understood,  and  are 
yielding  attractive-enough  designs,  to  warrant  serious  attention  from 
designers.  One  needs  only  to  view  tl;a  digital  filter  design  problem  as 
a problem  of  statistically  designing  a linear  minimum  mean-squared  error 
filter  or  predictor.  Elegant  solutions  abound.  An  added  dividend  paid 


by  this  way  of  thought  is  that  a class  of  digital  filter  designs  becomes 


4 


11.  Beginning  the  Design 

The  design  begins  with  specification  of  a bandllmlted  power  spectrum 
G(f)  of  the  type  illustrated  in  Figure  la.  It  may  be  helpful  to  visualize 
the  spectrum  on  the  passband  interval  -W  < f < W as  a signal  spectrum  S(f) 
and  the  spectrum  on  the  remainder  of  the  frequency  interval  (the  stopband) 
as  a noise  spectrum  N(f).  See  [5).  Then  S(f)  corresponds  to  the  desired 
passband  characteristic  and  N(f)  corresponds  to  the  desired  stopband 
characteristic.  Center  frequency  f^  (f^  = 0 in  Figure  1),  signal  band- 
width W,  stopband  rejection  n»  and  total  bandwidth  1/2T  are  design  param- 
eters. The  rejection  q may  be  visualized  as  a signal-to-noise  ratio 

n = A /A  . In  the  design  of  MA  filters  it  has  been  found  useful  to  include 
s n 

a transition  band  [5].  In  the  design  of  ARMA  filters  beginning  with  AR 
filters  we  have  not  discovered  a useful  way  to  specify  a transition  band. 

The  next  step  is  to  periodically  extend  the  spectrum  G(f)  using 
period  1/T  (or  foldover  frequency  1/2T)  and  scale  it  by  ^ to  obtain  a 
periodic  spectrum  G^(f): 


Gp(f) 


"iG(f)  , lf|  < 1/2T 


>^G(f  - m/T  - 1/2T  < f < m/T  + 1/2T  , m - 0,  ± 1,  ± 2,  . . . 


This  is  illustrated  in  Fig.  lb.  The  Fourier  coefficients 
1/2T 

S T / G (f)e~J^’'^‘'^df  , k = 0.  ± 1,  ± 2,  ...  (2 

-1/2T  P 

_ - 60  00 

specify  a covariance  sequence  with  = c The  sequence 

may  be  thought  of  as  the  covariance  sequence  of  a zero-mean  sampled 

00 

data  sequence  {x(kT)}  obtained  by  sampling  a continuous-time  random 


process  (x(t)}_^  that  has  bandlimited  power  spectral  density  G(f).  Then, 


I 


of  course,  = c(T=kT)  where 
1/2T 

c(t)  = / G(f)e~^^”^'^df  O: 

-1/2T 

00 

is  the  covariance  of  the  process  {x(t)} 

~00 

Another  connection  between  and  the  periodically  extended 

spectrum  G^Cf)  is 

G^(f)  - Ma)M(z-l)l,.,,p(^2.£T) 

M(z)M(z“^)  = E Cj^z’^ 
k=-«> 

00 

This  says  is  the  covariance  of  the  output  of  a digital  filter 

M(z)  that  is  driven  by  a white  sequence  with  unit  variance  and  that  the 
periodically-extended  spectrum  G^Cf)  is  the  magnitude-squared  frequency 
response  of  the  digital  filter  M(z). 

The  periodic  spectrum  G^(f)  Illustrated  in  Figure  1 is  a legitimate 
discrete-time  spectrum.  However,  it  is  not  rational  and  is,  therefore, 
not  the  spectrum  of  an  autoregressive  moving  average  digital  filter. 

This  is  another  way  of  saying  the  M(z)  of  (4)  is  not  a transfer  function 


for  a filter  that  can  be  realized  using  a finite  number  of  multipliers, 
adders,  and  delays  (memory).  Thus  the  approximation  problem  is  one  of 
designing  an  ABMA  filter  H(z)  whose  rational  spectrum  H(z=e^^’^^^) • 

H(z  =e^  ) approximates  the  irrational  spectrum  G^Cf).  Of  course, 

H(z=e'^ ^^^'^)H(z  is  simply  the  magnitude-squared  frequency 

response  of  the  filter  H(z).  This  brings  us  to  the  next  step  in  our 


design  procedure. 


III.  Designing  an  Autoregressive  Approximation 


The  problem  now  is  the  following:  given  a covariance  sequence  {c  }°° 

k 

corresponding  to  the  desired  spectrum  G (f),  find  a realizable  filter  H(z) 

P 

whose  periodic  spectrum  well-approximates  G (f).  Call  {r,  }°°  the  covari- 

p k 

ance  sequence  of  the  filter  H(z).  The  following  relationships  are  in 
force: 


H(z)H(z“^)  = Z r z 

k=-oo  ^ 

" 2^  H(z)H(z"^)z'*^"^dz  (5) 

j2TTfT.„,  -1  j27rfT.  “ -j2irfkT 

H(z=e-'  )H(z  =e-^  ) = ^ ® 

k=-co  ^ 

It  is  clear  that  for  purposes  of  matching  spectra,  the  covariances  come 
into  play  in  a more  fundamental  way  than  do  the  unit  pulse  sequences. 

In  (5)  the  contour  C lies  within  the  annulus  of  uniform  convergence  of 
H(z)H(z"^). 

One  approach  to  the  approximation  problem  is  to  design  an  AR  (N^) 
filter 

H(z) ^ (6) 

^ —0 
1-  Z a z 

1=1  ^ 

in  such  a way  that 


9 k Oyily.a.i  (7) 

N N 

For  large  values  of  N, , the  matching  of  {r,  and  {c,  }.  provides 

1 k 0 k 0 

effective  approximation  of  Gp(f)  by  H(z=e^^^^^)H(z~^=e^^^^^) . See 
equations  (2)  and  (5).  The  result  of  (7)  is  achieved  with  the  filter 
of  (6)  by  solving  the  normal  equations 
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Figure  2 shows  the  results  of  several  AR  designs.  Each  design 
was  begun  with  a G^Cf)  of  the  form  illustrated  in  Fig.  1.  The  idealized 
Gp(f)  is  superimposed  on  each  design.  Important  parameters  are  included 
in  the  figures. 

There  is  one  other  interpretation  of  the  H(z)  given  in  (6)  that  is 
worth  noting,  even  though  it  has  been  noted  many  places.  The  filter 


L(z)  = r a z~ 

1=1  * 

with  the  given  by  (8)  is  the  linear  MA  filter  of  order  N^^-l  that 

minimizes  the  mean  squared  prediction  error 


2 2 
C = E[x(mT)  - la.  x((m-il)T)]^ 

J!,=l 

00 

Here  {x(dT)l_^  is  the  sampled-data  sequence  that  has  spectrum  0^(1). 
The  filter 


K(z)  - 1 - £ a z 

1=1 

is  the  corresponding  "whitening"  filter.  Of  course  H(z)  = 1/K(z). 


Therefore,  to  the  extent  that  K(z)  whitens  G^Cf),  H(z)  inherits  the 
spectral  properties  of  G^Cf).  It  is  well-known  that  an  H(z)  designed 
as  outlined  above  is  stable  [9]. 
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IV.  Obtaining  Unit-Pulse  and  Covariance 
Sequences  from  the  AR  Design 

No  one  wants  a hlgh-order  AR  filter,  regardless  of  Its  frequency 
response  characteristics.  However,  the  hlgh-order  AR  design  serves 

another  very  useful  purpose:  It  provides  a ready-made  generator  for 

00  00 
a covariance  sequence  and  a causal  unit  pulse  sequence 

that  are  consistent  In  the  sense  that 


k = 0,  1,  2,  ... 


The  appropriate  equations  for  the  generation  of  the  unit  pulse 


00 

sequence  are 


k < 0 


k > 0 


Here  the  coefficients  are  the  feedback  coefficients  of  H(z).  See 

equations  (6)  - (8).  The  sjrmbol  6^^  denotes  the  Kronecker  delta. 

00 

The  generating  equations  for  are  not  much  more  difficult. 

The  following  linear  relationship  holds: 


So,  given  the  N covariances  {r.}  , for  example,  one  may  solve  for 

r^j  and  so  on.  Of  course,  the  are  available  as  r^  = c^, 

1 ■ 0,  1,  ...,  N^,  from  the  AR  design  procedure  of  Section  III.  For 
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completeness  we  show  In  Appendix  I how  to  solve  for  {r  }_  in  terms  only 
N 

of  the  {a,},  , and  thereby  completely  characterize  the  covariance  sequence 

JL 

00 

{r„}  . We  re-iterate  that  the  calculation  of  Appendix  II  is  not  required 

Ni 

in  our  design  procedure  because  the  are  available  as  r^^  = c^, 

1 = 0,  1,  . . . , with  the  Cj^  given  by  (2). 

Here  is  where  we  stand:  we  have  at  our  disposal  a systematic 

technique  for  generating  a high-order  AR  (N^)  filter  that  "approximates" 
the  idealized  spectrum  G (f).  This  AR  (N, ) filter  is  characterized  bv 
its  impulse  and  covariance  sequences  and  which  are  easily 

obtained  as  outlined  above.  The  sense  of  the  approximation  is  that 

— Cj^,  k — 0,  1,  2,  ...> 

In  the  following  section  we  use  the  algorithm  of  Mullis  and  Roberts 
to  approximate  this  generally  high-order  AR  (N^^)  filter  with  a much 
lower  order  ARMA  (N,M)  filter. 


V.  Designing  with  the  Mullis-Roberts  Algorithm 
Call  H(z)  the  transfer  function  of  our  AR  (N^)  design.  Let  H(z) 
be  the  transfer  function  of  an  ARMA  (N,M)  filter: 


H(z) 

Q(z) 

A(z) 


Q(z) 

A(z) 


^ ^ -M 

qQ  + qiZ  + . . . + q^z 


1 j.  -2  ^ ^ -N 

1 + a,z  + a_z  + ...  + a„z 
i z N 


(16) 


The  problem  studied  by  Mullis  and  Roberts  is  one  of  minimizing  the  I 

approximation  error  i 


-1/2T 


(17) 
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In  the  general  case  H(z)  is  arbitrary.  For  our  purposes  it  will  always 

2 

be  the  transfer  function  of  a high-order  AR  (N^^)  filter.  The  error  e 
may  be  written 

= T / ^ | ^df  (] 

-1/2T 

The  main  virtue  of  this  error  measure  is  that  it  leads  to  a tractable 
minimization  problem  [ 2 ] . It  is  a straight-forward  exercise  to  define 

A(z)  = H(z)A(z)  - Q(z),  invoke  Parseval's  Theorem,  and  use  (13)  and  (18) 
to  write 


2 

e = 


N N 


M N 


M 


1=0  m=0  “ I*'  k=0  i-0  ^ k-0  ^ 


(19) 


The  minimization  problem  is  one  of  minimizing  c with  respect 

M N 

to  the  ARMA  parameters  and  subject  to  the  constraint  that 

a^  = 1.  Write  this  as 


min  [e  - 2X(a'ljI  - 1)] 
a,q 

where  X is  a Lagrange  multiplier  and 


(20) 


® * A’  ®N-1’  ®0^ 


(21) 


ijj  = (0,  0,  . . • , 1) 

Note  a'iji  • a^. 

First  minimize  with  respect  to  q.  The  constraint  plays  no  role 
and  the  result  is 


N 


k = 0,  1,  . . . , M 


(22) 


i-0 


with  the  corresponding  to  H(z)  and  the  yet  to  be  determined. 
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Substitute  this  result  Into  (19)  to  obtain 
2 

= a Ra 

M 
I 

k=0 


e = a"Ka 
K = {r 


(23) 


*'-m|  . 


Here  the  notation  K = {d„  } denotes  that  d.  is  the  (£,m)-th  term  of  the 

£,m 

(N+1)  X (N+1)  matrix  K.  The  matrix  K is  positive  semi-definite  if  and  only  if 
{r,  }”  and  {h,  }„  are  consistent,  as  they  are  by  virtue  of  our  construc- 
tion  method. 

Now  our  constrained  minimization  problem  yields  the  following 
solution  for  a: 


Ka*  = xij; 


(24) 


Invoke  the  constraint  to  get 

a = (kjj/kQ,  kjj_j^/kQ,  1) 


(25) 


r-l 


where  k * (k^^,  . . . , k^)  is  the  last  column  of  K . Substitution  of 
a*  into  (22)  yields  q*  and  completes  the  design  of  the  ARMA  (N,M) 
filter 

'‘0  ■ '*1' 


H(z) 


+ q,Z  + ...  + q^ 


- , * -1  . 

1 + SqZ  + . 


^ * -N 

+ a^z 


(26) 


Call  {h,  }_  the  unit  pulse  sequence  and  {r,  } the  covariance 
k u k 


sequence  corresponding  to  H(z).  Since  a*  is  simply  the  solution  to  a 

normal  equation  Involving  positive-definite  K,  it  follows  that  {h,  } is 

k 

absolutely  summable.  That  is  H(z)  is  stable.  Further,  as  shown  by  Mullis 
and  Roberts,  hj^  = h^,  k * 0,  1,  . . . , M.  The  covariances  r^^  and  are  not 

All  of  the  results  of  this  section  are  contained  in  [2]. 


equal.' 
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VI.  Systematic  Design  Steps 

The  design  steps  are  summarized  below  and  Illustrated  In  Figure  3. 

1.  Specify  a periodic  power  spectrum  G (f). 

2.  Invert  G (f)  to  get  a design  covariance  sequence  (c.  . 

p IC  V 

3.  Solve  the  normal  equations  Ca  = c to  obtain  an  AR  (Nj^)  approxi- 
mation. Typically  choose  large  (e.g.  = 256). 

4.  Generate  ^ ’ generate  using  the  AR  gener- 

ating equation  (14). 

5.  Construct  the  matrix  K. 

6.  Solve  for  a*  and  q*. 

A 

7.  Scale  H(z)  to  get  desired  dc  response  Hq(z*1)  = H^. 

Actually,  steps  5 and  6 are  replaced  by  a Levinson-type  algorithm 
presented  by  Mullls  and  Roberts  for  efficiently  obtaining  a*  and  q*. 

The  FORTRAN  listing  of  Appendix  II  describes  & program  written  for  a 
CDC  6400  computer  with  SCOPE  compiler  to  implement  the  design  seeps  above. 
The  program  was  used  to  generate  all  of  the  designs  that  follow  In  Section 
VII.  The  number  of  storage  locations  required  for  the  TPLSV  subroutine 
to  solve  (8)  is  3Nj^.  The  number  of  multiplies  is  approximately  2N^  and  the 
number  of  adds  Is  approximately  N^ . The  number  of  storage  locations  re- 
quired to  solve  for  a*  using  the  Mullis-Roberts  algorithm  Is  0(N) . The 

2 

number  of  multiplies  Is  0(N  ).  The  notation  0(N)  Indicates  here  that 
0(N)/N  Is  a constant. 
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VII.  Designs 

The  designs  presented  in  Figures  4-8  all  have  a design  cutoff 
frequency  of  0.2Hz  and  a foldover  frequency  of  l.OHz.  This  does  not 
mean  the  3dB  point,  for  example,  is  under  careful  control.  It  only 
means  the  ideal  spectrum  with  which  the  design  was  begun  had  a cutoff 
frequency  of  0.2Hz.  Our  convention  in  the  figures  is  that  ARMA(N,M) 
always  denotes  a filter  designed  according  to  the  procedure  of  Section 
VI.  Of  course,  the  Chebyshev  and  Butterworth  design  of  Figures  7 and  8 
are  also  autoregressive  moving  average  filters.  These  filters  are  denoted 
CH(N,N)  and  BW(N,N),  respectively. 

The  results  of  Figure  4 illustrate  how  low-order  ARMA  filters  may  be 
used  to  approximate  high-order  AR  filters.  Figure  4(a)  illustrates  that 
an  ARMA  (8,8)  filter  is  a very  effective  approximation  to  an  AR(64)  filter 
when  only  30dB  stopband  rejection  is  desired.  The  complexity  of  the  ARMA 
(8,8)  is  one-fourth  the  complexity  of  the  AR(64)  and  the  passband  ripple 
of  the  ARMA  (8,8)  is  very  low.  Figure  4(b)  illustrates  that  an  ARMA  (32,32) 
filter  is  an  effective  approximation  to  an  AR(256)  when  90dB  rejection  is 
desired.  Figure  5(a)  Illustrates  what  happens  when  one  fixes  the  number  of 
poles  (N=8)  and  the  number  of  zeros  (M=4)  in  a relatively  low  complexity 
design  and  then  Increases  the  stopband  rejection.  At  30dB  rejection  the 
filter  characteristic  is  smooth  in  the  stopband  and  cuts  off  sharply.  At 
90dB  rejection  the  ripple  in  the  passband  is  severe.  Figure  5(b)  illus- 
trates that  if  additional  zeros  are  added  to  increase  the  complexity,  con- 
trol is  maintained  over  the  passband  ripple  even  as  the  stopband  rejection 
is  increased  to  90dB. 

Figure  6 shows  that  when  complexity  is  already  high  (N=16,  M^8)  and 
the  stopband  rejection  moderate  (60dB) , the  addition  of  zeros  does  not 
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dramatically  influence  the  magnitude-squared  response. 

Figure  7 shows  ARMA  (8,M)  approximationa  to  ideal  lowpass  spectra 
and  comparisons  with  Butterworth  and  Chebyshev  designs  of  order  8 and 
complexity  16.  If  only  30dB  stopband  rejection  is  desired  then  an  ARMA 
(8,4)  filter  can  be  designed  that  yields  smoother  passband  character- 
istics than  the  Cheybshev  filter  and  sharper  cutoff  than  either  the  Cheby- 
shev filter  or  the  Butterworth  filter.  Of  course,  the  passband  ripple 
is  larger  than  that  of  the  Butterworth  filter.  The  ARMA  (8,4)  filter  is 
less  complex  than  the  Chebyshev  and  Butterworth  filters  by  a factor  of 
3/4.  As  rejection  is  increased  to  60dB,  the  passband  ripple  of  an  ARMA 
(8,8)  Increases,  but  remains  smaller  than  the  Chebyshev  ripple  every- 
where in  the  passband  except  near  cutoff.  The  cutoff  is  comparable. 

Figure  8 shows  ARMA  (16, M)  approximations  to  ideal  lowpass  spectra 
and  comparisons  with  Butterworth  and  Chebyshev  designs  of  order  16  and 
complexity  32.  If  only  60dB  stopband  rejection  is  required,  then  an 
ARMA  (16,8)  provides  cutoff  comparable  to  that  of  the  Chebyshev  with 
smaller  ripple  except  near  cutoff.  For  90dB  rejection  an  ARMA  (16,16) 
provides  cutoff  comparable  to  the  Chebyshev  filter  of  the  same  complexity 
with  less  passband  ripple  except  near  cutoff. 

In  all  of  our  designs  run  to  date,  well- conceived  ARMA  designs 
exhibit  relatively  little  passband  ripple,  except  very  near  the  cutoff 
frequency,  and  cutoff  characteristics  comparable  to  Chebyshev  character- 
istics. When  filter  complexity  is  high  and  stopband  rejection  low, 
then  the  filter  characteristics  can  approximate  ideal  lowpass  character- 
istics much  more  closely  than  equivalent-complexity  Butterworth  or 
Chebyshev  designs.  Reference  the  ARMA  (8,8)  filter  of  Figure  4(a),  the 
ARMA  (8,12)  - 30dB  filter  of  Figure  5(a),  and  the  ARMA(16,24)  filter  of 
Figure  6. 
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VIII.  Conclusions 

The  results  presented  here  are  only  representative  of  what  one  may 
achieve.  Extensions  to  bandpass,  high-pass,  band  reject,  and  comb  fil- 
ters are  straightforward.  Several  points  should  be  made:  (i)  there  is 

no  requirement  in  these  designs  to  begin  with  a rational  analog  filter; 
(11)  there  Is  nothing  fundamental  about  the  shape  of  the  G(f)  we  have 
started  with  - an  inventive  designer  may  experiment  with  his  own  choices; 
(iii)  the  procedures  presented  here  permit  the  designer  to  exercise 
implicit  control  over  center  frequency,  bandwidth,  and  stopband  rejec- 
tion; and  (Iv)  one  may  systematically  design  filters  with  arbitrary 
numbers  of  poles  and  zeros. 

The  results  of  Section  VII  do  not  say  the  statistical  designs  of 
this  paper  are  better  or  worse  than  more  classical  designs.  Only  differ- 
ent. The  new  designs  do,  however,  offer  much  more  flexibility. 

We  speculate  that  the  fitting  of  long  AR  models  and  subsequent 
approximation  with  the  Mullls-Roberts  algorithm  may  provide  a useful 
method  of  fitting  ARMA  spectral  models  to  random  data.  The  question  of 
order  determination  would,  of  course,  be  crucial  here  as  It  Is  in  AR 
model-fitting  [6]. 


APPENDIX  I;  Generating  a Covariance  Sequence  from  an  AR  Model 


Write  out  (15)  for  k • 1,  2,  . . . , and  order  the  equations  in 
the  following  way: 

Qr  = r 


0 02  “3 


“1  °2  “3 


“2  “4  “5  \ ° ° 


“Ni  “Nj^-1 


“1  ° 


^ . (r^.  r^,  ....  r^,^) 


The  matrix  Q is  generated  as  follows:  begin  with  the  first  row 

(0,  a,,  0-,  ...  a„  );  left-shift  this  row  (n-1)  times  and  add  the 
J.  L 

mth  overflow  to  the  (n,m+l)  term  to  get  the  nth  row.  For  example, 
the  fourth  row  is  generated  in  the  following  way: 

0 a^  02  03  a^  05  . . . 

0 03  02  03  a^  03 

“3  “2-^4  “I'^S  “6  •••  “N3  ° 

The  vector  of  covariances  r is  obtained  by  solving  the  eigenvalue 
problem 

qI  - \\ 


for  the  eigenvector  I (»  r)  corresponding  to  A = 1 


<i«'sc'»oo  Of  .< '.ofiC  oo  f )or>oof  )Of  >«  70000000  r.  oooofjoooooo  or^fiocio 
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Appendix  II:  FORTRAN  Program  for  Design  of  ARMA  Filters 


PROGRAM  ARF  ILT<  INPUT,  OUTPUT  •FIL*«PL»T  ftPri.  = I r.PLT,  TAT  ri  =01JTPUT  ) 

PROGRAM  FOR  STATISTICAL  DESIGN  OP  AR  Afj[  ARMA  riGITAL  FILTERS 
WRITTEN  BY  J.  LUP.Y,  COLORADO  STATE  UMV.,  JINTER,  197F-77, 
PROGRAM  SET  UP  FO"  LOUPASS  E'ESKAS  AT  PEE  SENT  (AP  CR  AK^  A). 

DATA  READ  IN  SETS  OF  A rsROS.  VARIAPLL  ’jCC  TELLS  PRCGkAr-' 

THE  NUMBER  OF  G CARD  SETS  TO  READ  iN.crrt  pIFST  ap  STfTI'^ENT) 

INTEGER  PS 

INTEGER  NMRUO)  fNEM(in> 

INTEGER  TOOL(IO) 

DIMENSION  HQ(257) 

DIMENSION  E(257) 

DIMENSION  H<257),AC257),Q(?B7) 

DIMENSION  QN(257) 

DIMENSION  PN<257> ,RS<257) ,R  f 2B7) ,PP(257) 

DIMENSION  UAR(257) 

DIMENSION  PHAR<hl2),HAR(512>,WSa(512),OHAS»-(512),FRrC(512) 


ALL  DATA  CAROS  MUST  BE  PRESENT  EVEN  IF  THFY  ARE 
NOT  APPLICABLE  TO  THE  PPESErJT  RUN. 

BLANK  CARDS  MAY  BE  INSERTED  RQ"  N/A  LATA  CARCS. 

(FOR  AP  DESIGNS,  A BLANK  CAR?  IS  INSERTfO  FOR  T'JP 
NUMERATOR  ORJER  CARDCCARD  " IN  THE  SET  PF  G)) 

NOTE.,. EACH  DATA  SET  ALLOWS  DESIGN  OF  AT,  ARMA  0=  AR, 
BUT  NOT  BOTH. 

NOC=NUMBER  OF  DATA  SETS 
NCC  .LE.  10 
READ  (B,250)  NOC 

DO  2P0  IX  = l.NDC 


CARDl 


OPTION  CARD 

THIS  CARD  IS  READ  IN  1011  FORMAT 

A 1 OR  0 IS  PUNCHED  IN  C'^LUMNS  1-P,  WlTP  THE  PtLLCWI».r-  MLANlNG- 
1=DESIGN  AR  FILTER  0 = NO  AR  DC  SIGN 

1=PL0T  AP  PHASE  0=NO  RIOT 

1 = DUMP  RELEVrfiT  COERS.  P=NO  CU**? 

1=DESIGN  ARMA  FILTER  0=NO  ARUA  DESIGN 

IsPLOT  ARMA  PHASF  0=N0  PLOT 

1=DUMP  ARMA  COEF.  0=N0  DUMP 

1=PL0T  MAG.  SQUARED  0=NO  MAG.  SQ.  PLOT 

isPLOT  ON  MICROFILM  0=NO  PLOT 


ALL  COEFFICIENTS  IN  DU-ip  ARE  NORMALIZE?  FOR  0 D.C.  ON  pL'*T. 


READ  <5,270)  (T(;OL(I),I  s 1,10) 


CARD2 

JJAsNUMBLR  OF  DIFFERENT  FILTER  CRDERS  TC  PE  OESICNE? 
UNDER  THE  ACTIVE  OPTION  CARD, 

JJA.LE.9 

NOTE  THAT  JJA  IS  THF  OP  lOOF  ?Aca*'ETfe  I'  1‘^E  R ^ A p 
STATEMENTS  OF  CAROS  H AN'I  (, 

FORMAT  111 
READ  (5,2b0)  DJA 


CAR03 

bWS=SlGNAL  EmMUWIDTH  IN  H7  3UN  = 'i(I£R  PW  I"  H7 

dWS^bwN  MUST  PF  .LE.  FCLIOVER  fr^o, 

FORMAT  3F3.2 
READ  (5,280)  BWS.pwn 


'Jd 


OOOOO  DD  oonooooooo  TOO  o oooooo  ooooo  ooooooo 
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CARO<» 

lONE  IS  THE  SMALLEST  S/M  RATIO.  ITVt  IS  THE  LARGEST  S/N. 
ITHREE  IS  THE  S/N  LOOP  INCREMENT. 

A TYPICAL  CARD  A IS  AS  FOLLOWS. .. 0030 0 90 03 ... THI S WOULD 
CAUSE  FILTERS  WITH  S/N  RATIOS  OF  30*60*90  DECIBELS  TO 
BE  DESIGNED. 

FORMAT  313 

READ  (5*260)  lONE *I TWO* ITHREE 
CAR05 

NMR  ARRAY  TAKES  ON  THE  VALUES  OF  NUM.  ORDER  FOR  ARMA  FILTERS. 
EACH  VALUE  OF  NMR  SHOULD  BE  .LE.  128. 

FORMAT  1013 

READ  (5*350)  (NMR(I)*I  = 1*JJA) 

CARD6 

NEN  ARRAY  STORES  VALUES  OF  DEN.  ORDER  FOR  AR  OR  ARMA. 

A MAX.  OF  10  VALUES  CAN  PE  STORED.  EACH  OF  THESE  MUST  BE 
.LE.256  FOR  AR*  AND  .LE.126  FOR  ARMA. 

FORMAT  1013 

READ  (5*350)  (NEN(I)*I  - 1*JJA) 

PI  = 3.141592659 
T 5 0.5 
PPP=0.5/T 

T IS  THE  SAMPLING  PERIOD. 

PPP  IS  THE  FOLDOVER  FREQUENCY. 

MAUG  = 1024 
NN  = 512 

THE  USER  SHOULD  DEFINE  VARIABLES  APPLICABLE  TO 
THE  TYPE  OF  FILTER  TO  BE  DESIGNED. 

THE  VARIABLES  FO*PWS *BWN*FQ . . . DEFINE  THE  IDEAL 
POWER  SPECTRAL  DENSITY  OF  A LOWPASS  FILTER. 

THE  FOLLOWING  SECTION  BRACKETED  BY  STARS  NEEDS  TO  BE 
CHANGED  FOR  OTHER  TYPES  OF  FILTERS. 

THIS  SECTION  COMPUTES  THE  COVARIANCE  OF  THE  IDEAL  P.S.D. 


BWQ=PPP>(BWN-»BWS) 

FO=PPP-BWN/2.0 

FQ=(PPP-BWN^BWS)/2.0 

COMPUTE  COVARIANCE -RS=SIGNAL  CDV.*  PN=NCISE  COV. 
QN=TRANSITION  BAND  COVARIANCE 
PN(1)  = 2.0  • BWN 
RS(1)  = 2.0  * BWS 
QN(1)  = 0.0 
DO  100  J = 2*257 

RS(  J)=SIN(2.*PI*FL0AT(J-1)*T*RWS)/(PI*T*FL0AT<  J-D) 

PN< J)=2.0*SIN(P1*BWN*FLOATC J-1) ♦T>/(Pl*FLOATt J-1)*T)*COS(2.0* 
CPI*FO*FLOAT(J-l)*T) 

QN(J)  = 0.0 
100  CONTINUE 


AN^NOISE  LEVEL.  AS=SISNAL  LEVEL(A  VARIABLE) 

aqstrans.  band  level 

AN  = 1.0 

DO  230  KK  = IONF*ITWO*ITHREE 
AS  = 10  * * KK 

DO  230  NNA  : l.JJA 
NOEN  : NEN(NNA) 

MAR  s NOEN  ♦ 1 
IF  (TOOL(l).EQ.l)  60  TO  110 
C FOR  ARMA  DESIGNS*  AR  0R0ER=?56 

MAR  = 257 
MR  = MAR  - 1 

R(l)  s AS  * RS(l)  ♦AN  * PN(l)  *40  • QN(1) 

DO  120  J = 2.MAR 

R(J>  s AS  * RS(J)  ♦AN  * PN(J)  ♦ AQ  * QN(J) 


110 


ooonnoo  o o r)r>  nrx^ofjo 


R=TOTAL  COVARIANCE  VECTOR.  EQUALS  SIGN^L+NCISE  COV..  .EISHrED 
BY  AS  AND  AN. 

ARRAY  3N  ALLOWS  THE  DESIGNER  TO  SPECIFY  A TRANSITION 
BAND  COVARIANCE  WITH  WEIGHT  AQ  AND  CENTER  FREQ.  FG. 

120  CONTINUE 

BS  = MAK  - 1 
DO  130  LQ  = 1*BS 
LOl  = LQ  ♦ 1 
PPCLQ)  = - R<LR1) 

130  CONTINUE 

CALL  TPSLV  (BS»R*PP»WAR,E) 

SUBROUTINE  TPSLV  SOLVES  T0FPLIT7  SET  OF  EOUATIONS- 
R*WAR=PP. 

JS  = MAR 
lAO  JT  = JS  - 1 

WAR( JS)=WAR( JT) 

IF  <JS.EQ.2)  GO  TO  150 
JS=JS-1 
GO  TO  IRO 
150  CONTINUE 

DO  160  I = 1<?57 
A( 1 ) = O.C 
Q(I > = 0.0 
160  CONTINUE 

WAR (1)  =1.0 

GSQR=AR  FILTER  GAIN  COEF.  SQUARED 
GSOR  = R(l) 

DO  170  JP  = 2,I^AP 
6SQR  = GSQR  ♦ UAR(JP)  * R<JP) 

170  continue 

Q(l)  = GSQR  * * .5 

GENERATE  IMPULSE  RES®.  SEQUENCE  FOR  MULl  I S/R  CPEP  TS  ALGOPITH'I. 
CALL  IMPT  <257,H,WAR*Q> 

IF  (TOOLU)  .NE.l)  GO  TO  180 

THL  FOLLOWING  STARRED  SFCTI'-N  IS  USED  p(R  AR  DESIGNS  ONLY. 
ALTHOUGH  A R < V ) =A RMA ( N , 0 > ♦ A TIME  SAVINGS  RESULTS 
BY  HAVING  A SEPARATE  AR  DESIGN  ’’ROCriURf.  T^E  AR 
SOFTWARE  IS  USED  FOR  AR‘<A  DESIGNS  ALSO. 

*•*«*****•***♦***•************«***••*****•*«*****«***«*«* 

CALL  ARMAMG  ( Q , WAR ♦MAR ♦ WSQ »PH ASE  , 1 ) 

CALL  SPLPLT  ( NN  ♦ WSQ  tPHASE  ♦»<  AR  ,PH AR  tFRE Q ) 

WRITE  (6t330> 

WRITE  (6^310>  NOFN^kk 
WRITE  (6f330) 

WRITE  (6«3R0)  BUSfBWN 
WRITE  (6f330) 

IF(TOOL<7).EQ.O)  GO  TO  175 
KKL  s TOOLI2) 

CALL  PLOT  <HARfPHAR^NN$KKLtTOOLfFREQ*PPP) 

175  CONTINUE 

WRITE  (6f33D) 

IF  <T0DL(3)  .EQ.O)  GO  TO  IPO 
WRITE  (6$290) 

WRITE  (6f330) 

CALL  PRINT(WARfR^H^ AfMARt3) 

C 

0 «*•******•**«*.*•*«*•••«•***«•*•****•*»«•*****.*.*•*.•*** 

180  CONTINUE 

C GO  TO  240  IF  AN  AP  FILTER  T**  3EING  fESIGNEP. 

IF  (TOOL(4 ) .FO.' 1 SO  TO  230 
Q(l)  = 0.0 
DO  190  MY  = 

HQ(MY)  = ^^(MY) 

190  CONTINUE 

IDIF  = NEN(NNA)  - NMR(NNA) 

JDIF  = 0 

IF  (IDIF.LT.O)  JDIF  r 1 
IF  (IDIF.LT.O)  IDIF  = - 101^ 

IF  (IDIF. EQ.O)  GO  TC  2P0 


3 


oooo  o ooo 
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GENERATE  ARMA  COEFFICIENTS  GIVEN  IMP,  AND  COV,  SEQUENCES. 

CALL  MULRB  <MAR*H.R,A) 

CALL  NUMCF  ( Ot LSTAR • MAR ♦HQ ♦ A ) 

CALL  IMPT  (257,H,AtQ) 

WRITE  (6«330) 

CALL  ARMAMG  < Q ♦ A *M AR lUSQt PH ASE ♦ LST AR ) 

CALL  SPLPLT  < NN^USQ ♦PH ASE *H AR ♦ PHAR ♦ F REQ) 

WRITE  fSt330) 

WRITE  (G«330) 

IF  CIOIF.NE.O)  GO  TO  210 
NUMORD  s NFEN 
60  TO  220 
210  CONTINUE 

NUMORD  = LSTAR  - 1 
220  CONTINUE 

WRITE  (Gt320)  NFEN^NUMORD^KK 
WRITE  (6t330) 

WRITE  (6f3A0)  BUSfBWN 
WRITE  (6^330) 

KKK  = T00L(5) 

IF  (TOOL! 7).E0.0)  GO  TO  225 

CALL  PLOT  (HARfPHARfNNtKKKfTOOLfFREQfPFP) 

225  CONTINUE 

WRITE  (G^330) 

IF  (T00L<6)  .EQ.O)  GO  TO  230 
WRITE  (6f300) 

WRITE  (Gf330) 

IF  (MAR. LT. LSTAR)  MAR=LSTAR 
CALL  PRINT  (Q^AfR^HtMAR^ A) 

230  CONTINUE 
2A0  CONTINUE 
250  FORMAT  (111) 

260  FORMAT  (313) 

270  FORMAT  (1011) 

280  FORMAT  (3F3.2) 

290  FORMAT  (9Xfl3H  FILTER  COEF.,5XtlOH  COV.  SEC.,1?X^10H  IMP.  SEQ.) 

300  FORMAT  (9XtllH  NUM.  COEF . ♦ 1 1 X ♦ 1 1 H OEN.  COEF,^9XtlOH  COV.  SEQ.fllX. 
IlOH  IMP.  SEQ.) 

310  FORMAT  (13H  AR  FILTER  S/N  = 10**.13) 

320  FORMAT  (6H  ARM A ( ♦ I 3 ♦ IH ♦ ♦ I 3, H ) , 5X ♦?H  S/N  = lC**fI3) 

330  FORMAT  (IHO) 

340  FORMAT  (5H  BWS=fF5.2f5H  BWN=.F5.2) 

350  FORMAT  (1013) 

END 

SUBROUTINE  MULRB (N ♦ Hf R ♦ A) 

COMPUTES  DFN.  COEFS.  FOp  AR'-A. 

A ARRAY  CONTAINS  THESE  COPF. 

DIMENSION  A(N)fE(257).C(257> 

DIMENSION  H(N)fR(N) 

ALPH  = R(l)  - H(l>  * * 2 

A(l)  = 1.0 

B(l)  s 1.0/ALPH 

C(l)  = H(1)/ALPH 

OELT  = 1.0  ♦ H(l)  * * 2/ALPH 

NSTAR  = N - 1 

DO  100  I = 2#257 

B(I)  : 0.0 

C(l)  = 0.0 


100  CONTINUE 


iOtJOOc  > 


DO  170  I = 1*NSTAR 
GMA  = 0*0 
BET  = 0.0 
DO  110  K = 1«I 

BET  = BET  ♦ A(K)  * R(I  ♦ 2 - K> 

GPA  = GHA  - A(K)  * H<I  ♦ 2 - K) 

110  CONTINUE 

THET  = BET  * B(l>  * C(l> 

PHI  = BET  * C(l)  ♦ GKA  * nrUT 
ALPH  = ALPH  - BET  • THET  - GMA  * PHI 
IF  (ALPH)  130»l?0*i:^0 
120  WRITE  (6,190) 

CO  TO  180 

13C  DELT  = OELT  ♦ (PHI  * * 2)/aLPH 

IPLUS  =1*1 
DO  lAO  < = ItIFLUS 

A(K)  = A(K)  » BET  * B(1  ♦ 2 - K>  - GMA  * C(T  ♦ 2 - K) 
lAO  CONTINUE 

DO  150  K = 1, IPLUS 

E(K)  = B(K)  - (THET/ALPH)  * Ad  ♦ 2 - K) 

150  CONTINUE 

DO  160  K = 1. IPLUS 

C(K)  = C(K)  - (PHI/ALPH)  * A(I  ♦2-0 
160  CONTINUE 

A(1  42)=  0.0 
B(1  ♦ 2)  = 0.0 
C(I  ♦ 2)  = 0.0 
170  CONTINUE 
180  RLTURN 

190  FORMAT (3X.15H  ERROR  IN  MULR3) 

END 

SUBROUTINE  AR M AMG  ( 0 ♦ A ,M , WS9R  ♦ = •“  t LST  A R ) 

COMPUTES  ARMA  '"AG.  AND  PHASE  GIVEN  NUM.  AND  DEN.  CCEF  . 
THE  a ARRAY  CONTAINS  NUM.  GOEFS. 

THE  A ARRAY  CONTAINS  DEN.  COEFS. 

LSTAR=NUM.  ORDER  ♦! 

N=DEN.  ORDER  ♦! 

DP^ENSION  PH(512)  »WS.RR(512)  .WSPl  (512)  ♦ WS02(512).G(N)  , A(  N) 
COMPLEX  DATA1(102A).DATA2(1024) 

SUMl  = 0. 

SUM2  = 0. 

DO  100  I r 1,N 
SUM2  = SUM2  ♦ A(I) 

100  CONTINUE 

DC  no  1 = l.LSTAR 
SUMl  = SUMl  ♦ Qd) 
no  CONTINUE 

DO  120  1 = 1«N 
Ad)  = A(I)/SUM? 

120  CONTINUE 

DO  130  1 = l.LSTAR 
Qd)  = Jdl'SUMl 
136  CONTINUE 

DO  140  I = l.N 
DATA1(I)=CMPLX(A(I).0.0) 

140  CONTINUE 

DO  150  1 = l.LSTAR 
DATA2(1)=CMPLX(0(I),0.0) 

150  CONTINUE 
MP  = N ♦ 1 
DO  160  1 = MP,1024 
DATAKI  )=CMPLX(  0..0.) 

160  CONTINUE 

MP  = LSTAR  ♦ 1 

DO  170  I = MP,1024 

OATA2(1)=CMPLX(0..0.) 


170  CONTINUE 


onr>  oooo  nriooon 
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CALL  FOURT  ( D AT  A 1 , 1 024 ,1 , - 1, 0,0.0 
CALL  FOURT  (DAT A2 , 1 024 . 1 ♦ - T.OtO.O 
00  180  I = 1,512 

PH(I)=57.2957795*((ATAT42(AIt'AGOATA2(I)),REAL(rATA2(I)))>- 
CATAN2(A1HA6(0ATA1  (I)>,RFALnATAl(I>)>) 

PH(I)  = - PH(I) 

WSQKl)  = REAL(DATA1(  I ) * (OATAl  ( I > ) ) 

WSQ2(I)  s REAL(0ATA2(  1>  * CON JG ( D AT A 2 ( I > ) ) 

180  CONTINUE 

00  200  I = 1,512 
IF  (USQKI)}  190,200,190 
190  USQR(I)  s USQ2(I)/USQ1(T) 

200  CONTINUE 
RETURN 
END 

SUBROUTINE  SPEPLT (NN , USD ,PHA sr , H AR  , = H A R , FR LO ) 

COMPUTES  A SEQUENCE  OF  MAG.  AND  »HASE  PCINTS  TO  bE  SENT  TC  ’LOT. 
MAG.  SQUARED  VECTOR  IS  WSO. 

HAR  VECTOR  CONTAINS  WSQ  IM  flEClBELS. 

PHAR  CONTAINS  A SEQUENCr  OF  PHASE  POINTS, 

DIMENSION  WSQ(NN),PHASE(rjN)  ,HAP  ( NN  ) , PH  AP  ( S :< ) , Tr  f q { ';n  ) 

DO  100  I = 1,NN 

HAR(I)  = 10.  * ALOG10(USQ(I) ) 

PHAR(l)  = - PHASECI) 

FREQCIl  = FLOATCI  - D/512 
100  CONTINUE 
RETURN 
END 

SUBROUTINE  PR  TNT ( A , B , C ,0 ,NA .N®  ) 

PRINTS  OUT  ANY  ONE ,TWO, THR EE .OR  FOUR  1-riMENSICMAL, EQUAL 

LENGTH  ARRAYS.  NA  IS  THE  NUMBER  OF  ARRAY  ELEMENTS  TO  BE  PRINTEC. 

DIMENSION  0UT(4,257) 

DIMENSION  A (257), P( 257) ,C ( 257 ) , 0 ( 25 7 ) 

DO  100  J = 1,4 
DO  100  1 = 1,257 
OUT(J,I)  = 0.0 
100  CONTINUE 

DO  140  I r 1,NA 
IF  (NB.EQ.l)  GO  TO  130 
IF  (NB.EQ.2)  GO  TO  120 
IF  (NtJ.EQ.3)  GO  TO  110 
OUT(4,I>  = D(I) 

110  0UT(3,I)  = C(I) 

120  OUT(2,I)  s B(I) 

130  OUT(l,I)  = A(I) 

140  CONTINUE 

DO  150  1 = ],NA 

WRITE  (6,160)  I,(OUTf J, I),J  r 1,4) 

150  CONTINUE 
RETURN 

160  FORMAT  ( 2X , I 3, 5X , 1P4E2 0 . 7) 

END 

SUBROUTINE  TPSL V ( N ,R , PS, W, E ) 

ROUTINE  WRITTEN  BY  DAVE  FARDEN,  COLORADO  STATE  UNIV. 

DIMENSION  R (N),PS(N),W(N),r(N) 

NMl  = N - 1 
CINV  = 1.0/R(1) 

DO  100  1 = 1,NM1 
R(I)  = CINV  * R(I  ♦ 1) 

100  PS(  I)  = CINV  * PS(1) 

PS(N)  = CINV  * PS(N) 

W(l>  = PS(l) 

Ed)  = - R(l) 

AMB  = 1.0  * R(l)  * R(l) 


u w 
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DO  IfcC  I = 1*NM1 
II  = I ♦ 1 
THLT  r PS(Il) 

E.T  = - Rfll) 

DO  110  L = 1»1 
LI  = 1 - L ♦ 1 
THET  = THET  - U(L)  * RCLI) 
110  ET  = ET  - R(L)  * E(L1) 

Cl  = THET/AMB 
C2  = ET/APB 
11  = I ♦ 1 
IBAR  - 11/2 
DO  120  L = 1*1 
LI  = 11  - L 


120 


130 

lAO 

150 

1£0 


170 


* E(LI) 


E(LI) 
C3 


U<L)  = U(L)  ♦ Cl 
DO  130  L = 1,IBAR 
LI  = II  - L 
C3  = E<L) 

E(L)  = EtL)  ♦ C?  * 

E(L1)  = E(L1)  ♦ CP 
IF  (IBAR  - 1/2)  150,150»14C 
E(IBAR)  = C3  * (1,0  ♦ r?) 

U(I1)  = Cl 
E(Il)  = C2 
AK.B  = AMB  - CP  * ET 
UNNORHALIZE  R ADD  PS 

vJ  = N 

C = 1.0/ClDV 
DO  170  1 = l.NPl 
R(J)  = C * R(J  - 1) 

PS(I)  = C * PS(1) 

J = J - 1 
PS(N)  = C 
R(l)  = C 
RETURN 
END 

SUBROUTINE  PLOT ( PAR  ,PHAR ,NM ,T OOL »F REQ» PPP ) 

MAPA  AND  PAPN  ARE  PLOTTING  ROUTINES  AVAILABLE  AT  COL. 
MAPA  IS  A PAPER  PLOT. 

MAPM  IS  A MICROFILM  PLOT. 

INTEGER  TOOL(IO) 

DIMENSION  HAP(KN),PHAR(NN) 

DIMENSION  PREQ(NM) 

INTEGER  XT(B)fYT(8)*PT(P),^’T(7) 

XTd)  = lOHFREOtH?) 
lOH 

lOHMAGNSO (D&> 
lOH 

lOHPHASF  rPGR) 
lOH 


* PSCN) 


XT(  2) 
YT(1) 
YT(2) 
PTC  1) 
PT(2) 


DO  100  I = 3»B 
XT(1)  = lOH 
YT<I)  - lOH 
PT(I)  = ION 
100  CONTINUE 

MT(1)  = lOHFILTER 
MT(2)  = lOHOESIGN 
MT(3)  = lOH 
MT(4)  = lOH 
MT(5)  = lOH 
MT(6)  = lOH 
MT(7)  = Art 

CALL  MAPA  ( IfFRLQfHAR ,1 «NN.n, O.PPP, 
CALL  MA»A  (2,FREQ,HAR,1,NN,0.0,'’P3, 
CALL  MADA  (A*FRrOfHAR,l,NfSO.O.PPP* 
IF  CTOOL(8>.EQ.O)  GO  TO  110 
CALL  MAPM  <1,FREQ»HAR,1«NN, O.O^PPP* 
CALL  MAPM  (2*FRCO*HAR,l,Nf'»O.C.PPP* 
CALL  MAPM  (A«PRLO*HAR«1.NV.O.O.PPP, 


1P0.C.30.C.>T,YT,''T* 
120. 0*30. O^XT.YT^'^T. 
1?0.C*30.0*XT»YT,*-'T, 

lP0.r»30.0,XTfYT»wT, 
120. 0 *30. OfXTfYTtWT* 
1?C.0*30.0*XT.YT,vt, 


TATE  J. 


- 1) 
- 1) 
- 1) 

- 1) 
- 1) 
- 1) 
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110  IF  iN.EQ.O)  60  TO  1?0 

CALL  MAPA<l,FREQ,PHAR,l,NN*0.0*  = Po»-3f  r.».UO.«'<T,YT,'»T,-l> 

CALL  MAPA(2*FREO*PHAR*l,NN»0.0»PPP«-360,,3f 0.*XT,YT,MTf-l) 

CALL  MAPA<AtFREO*PHAR*l*NN*P.n*PpP*-360.*3<’  0.*XT*  YT*MT»-n 
IF  (TOOLf8).EQ.O)  60  TO  120 

CALL  HAPM(l»FREO*PHAR*l,NN*0.0fPpP,-360.»360.*yT,YT,MT,-l) 

CALL  MAPH(2*FREO*PHARtl,NN*0.0,PPP«-360.,3f>0.*XT,YT,PT,-l ) 

CALL  MAPHCR»FRF.Q*PHAR  tl,NN*0.0*PPP»-360.  P.,XT,YT*w|T*-l ) 

120  RETURN 
END 

SUBROUTINE  IMPT <N »Hf A • Q) 

C 

C CONFUTES  ARMA  IMP.  RESPONSE  GIVEN  COEFFICIENTS. 

C N VALUES  OF  THE  RESPONSE  ARE  COMPUTED  A> D STORED  IN  H. 

C N.LE.257 

DIMENSION  HCN)*A(N),Q(N) 

DO  130  1 = 1«N 
HU)  = 0.0 
SUM  = 0.0 
DO  110  J = 2«N 
IF  (1  - J ♦ 1)  120,120*100 
100  SUM  = SUM  - A<J)  * HU  - J ♦ 1) 

110  CONTINUE 
120  CONTINUE 

HU)  = QU)  ♦ SUM 
130  CONTINUE 
RETURN 
END 

SUBROUTINE  SHIFF (H,K ,NL,NX > 

C 

C THIS  ROUTINE  PERFORMS  A pra'^'^IAN  SHI'^T  CN  THE  I**PULf~ 

C SEQUENCE  H.  K IS  T^E  NU'^PER  OF  SHIFTS  TO  BE 

C NL  IS  THE  LENGTH  OF  THE  H SE3UENCE  SNP  IT  MUST  '^E  CCNSISTiNT 

C WITH  THE  LENGTH  OF  THE  H SEO.  SENERATFO  BY  IMPT. 

C NX  =1  FOR  LEFT  SHIFTS  AND  0 FOR  RIGHT  SHIFTS. 

DIMENSION  HCNL) 

IF  (NX.EQ.l)  GO  TO  130 
I = NL 

100  H<I)  = HU  - K) 

1 = 1-1 

IF-  U.GT.K)  GO  TO  100 
1 = 1 

IF  (K)  120«120«110 

^ 110  IF  U.GT.K)  60  TO  120 

I HU)  = 0.0 

^ I = I ♦ 1 

f GO  TO  110 

i 120  CONTINUE 

1 60  TO  150 

1 130  JX  = 257  - K 

n DO  140  I = 1,JX 

k!  HU)  = HU  ♦ K) 


i 


c 

c 

c 

c 


140  CONTINUE 
150  CONTINUE 
RETURN 
END 

SUBROUTINE  NUMCF «Q*LSTAR *MAR , HO*A) 

GENERATES  ARMA  NUM.  COFFS.  GIVEN  DEN.  CDCFS.  AND  IMP. 
SEQ. (UNSHIFTED). 

DIMENSION  Q(LSTAR)«HQ(MAR)«A(mar) 

NSTAR  = MAR 
DO  120  I = IfLSTAR 
DO  110  J = l.NSTAR 
IF  U - J ♦ 1)  120«120U00 
100  QU)  = QU)  * A(J)  * HQU  - J * 1) 

110  CONTINUE 
120  CONTINUE 
RETURN 
END 


nr>n 


t 
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100 

110 
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SUBROUTINE  COV < HQ, R,M  AR  * 1 D! «" ) 

MODIFIES  COV.  SEQ.  FOR  a L'‘FT  SHIFTED  IKP. 

DIMENSION  R(rS7),MQ(257) 

00  110  1 = IfMAR 

SUM  = 0.0 

00  100  J = IflOTF 

SUM  = SUM  - HQ(J)  * HOn  ♦ J - 1) 

CONTINUE 

R(l)  = R(I)  ♦ SUM 

CONTINUE 

RETURN 

END 


Sf  c. 
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